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A program generator to manipulate automatically Poisson series over tlie field of rational numbers is applied 
to develop the limit cycle of Van der Pol's equation in the powers of the small parameter. The results indicate thai 
the recurrence relations in what Melvin calls the algorithm of the shifted phase are stable. 
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In a recent note [l] 1 , Melvin presented a new set of rules to expand the limit cycle of Van der Pol's 
equation 

y + y = ey(l - y 2 ) 

in power series of the small parameter e. He carried out the calculations in single precision floating point 
arithmetic; he suggested thai they be repeated in exact form, i.e. that the coefficients in the series be 
produced as quotients of relatively prime integers. To achieve this result we applied to this problem a program 
generator of our own, called MAO, to manipulate rational Poisson series by computer. In the present note we 
show how this processor enables us to simplify drastically Melvins algorithm; then we present the actual 
expansions for the limit cycle and its frequency to e 8 . 

Melvin's rules are expressed in terms of a new independent variable 

T = col 

by means of which the original equation is transformed into 

0) 2 y" + y = ecoy'(l — y 2 ). 

The cycle and its frequency are expanded in formal power series 

7 = 2 -,yn(r)e n , (1) 

n>0 nl 

lo = J, -7*J n € w (2) 

n>0 nl 

of the small parameter e. Melvin showed that the following parity rules are consistent: 
(i) For n even, the term y n is an even function of r ; 

(ii) For // odd, the coefficient co n is and the term y n is an odd function of r. 
The terms y n and the coefficients a> n are to be determined recursively, starting with 

j = 2 cost and co = 1. 



* An invited paper. 
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1 Figures iu brackets indicate literature references at the end of this paper. 
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Let 

F = ea>y'{\ - y 2 ) - u 2 y' 
be expanded as the power series 



F = 2 ~ } e n F n 

n>o nl 



As we introduce the intermediate quantities 



a = j 2 , i' = toy ' 



n>o w i n>o w I 



we observe that the coefficients in the power series 

U = Z, — U n € n , v = 2j ~ v n* n i w = 2j ~ w n* n 

are given by the expressions 

2 ( ) ymJn-m, "w = 2 (_ 

0<m<w \ m / 
where ( J stands for the binomial coefficient n\/m\{n — m)\. Therefore 



, ^mj n—m 
0<ra<n \ m / 0<?n<w \mj 



T] V n -i ~ Zj \ VmUn-\-m ~ Zj ( I 

0<m<M-l \ m J J 0<m<n \ m / 



™mjn- 



We denote by F% the value which F n takes when we set y n = 0. In view of the parity rules, F n is an even or 
odd function of r respectively when n is even or odd. Now it turns out that the term y n in the limit cycle is a 
solution of the non-homogeneous linear equation 

7n + Jn = F* (3) 

Hence 

Jn = Jn + J* 

with the homogeneous general solution 

Jn = c n cos T + s n s i n T 

and y * a particular solution of eq (3). In application of the parity rules, we set s n — when n is even and c n 
= when n is odd. For the particular solution y* to be periodic in t, the right hand member F% must not 
contain terms in sin r and cos r. These critical terms will be removed from F* by following the periodicity 
rules: 

(i) When n is even, F * contains a term in cos r whose coefficient is of the form a — bco n . Therefore, 

we continue determining the frequency in the limit cycle by setting 0) n = a lb. The coefficient c n in 

the term y° n will be determined at the next order, 
(ii) When n is odd, F% contains a term in sin r whose coefficient is of the form a — bc n - 1 . The 

determination ofy n _ x is completed by setting c n _ x — a\b. We also sets n = which implies that y° n 

= in all components of the limit cycle at odd orders. 
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For our processor of Poisson scries, the terms y n are represented as pointers to list structures in which the 
trigonometric components are chained dynamically. The parity rules and die periodicity rules are all we need 
to know about the structures of the terms v„ in order to construe! the limit cycle. But, in Melvin's 
implementation, the Poisson series are represented by the static arrays of their numerical coefficients. Thus 
one needs to know that the term v, , is a Fourier sum of the form 



yn 



2j c ntTn cos(2m + 1 )r 



when n is even, and of the form 



Yn = X s n,m sin(2m + 1)t 



when n is odd. Expressed in terms of the matrices c n%m and s nm . the algorithm of Lindstedt-Poincare gives 
rise to a large set of intricate algebraic formulas. 

The tables on the next page present the development of the limit cycle to order 8. We confirm that the 
corresponding coefficients produced by Melvin in floating point format arc exact to 6 digits with rounding for 
the last digit. As the tables show, numerators and denominators in the coefficients of the cycle and of its 
frequency are rapidly increasing in absolute value. Our programs written in PL/ 1 for the IBM Optimizing 
Compiler (version 1. release 3.0) make use of the arithmetic hardware to operate on integers in decimal 
notation with at most 15 digits. Hence, the program encountered a fixed point overflow in the course of the 
operations at the 9-th order. A preliminary version in which the normalizing factorials were omitted in ( 1 ) and 
(2) raised the fixed overflow condition at order 5. 

Another way of determining the limit cycle uses complex coordinate's (//. / ). They are introduced by 

// = y — iy 
v = y + iy 
They transform van der PoFs equation into 

il=iu - l(u-r)((u+ v) 2 -4) (4) 

o 

A conjugate complex equation holds for?;. As the cycle in real form is a Fourier series with period 27t/(l> we 
introduce a new independent variable by 

£ = e i<ot 

The solution is then found as a formal power series of the form 

u = £ e"»n(0 

where the u r ,(^) are odd polynomials in ^. The frequency CO will be found at the same time in the form of the 
follow ing scries 



= 2 oj„e" 



The new independent variable gives rise to the differential operator 

^ d Id 

D = %—= 

dJ; ico dt 
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Table I. Frequency. 



n 


<*> n 





1 


2 


1 
8 


4 


17 
128 


6 


175 
6144 


8 


4752293 
884736 



Table II. Terms of even order in the cycle. 



>'4 



Te 



COS T 

cos 3r 
cos 5r 
cos It 
cos 9r 
cos 11 r 
cos 13 r 
cos 15r 
cos 17r 



1 

32 
3 
16 
5_ 
48 



23 
2048 
101 
512 
1865 
4608 
1379 
4608 
183 
2560 



2.58095 

1 1 79648 

120305 

196608 

1644175 

1769472 

10923199 

4423680 

1769369 

819200 

409871 

460800 

715247 

5160960 



1958726609 

141557760 

3919930525 

2.54803968 

147344494843 

31850496000 

161113663733 

5898240000 
1359229760383 

46448640000 
2076538440769 
130056192000 
526426361 



115605504 
392636471 
743178240 



Table III. Terms of odd order in the 



73 



I 



sin 3t 
sin 5r 
sin It 
sin 9r 
sin llr 
sin 13r 
sin 15r 



45 

256 

85 

384 

_7 

96 



3895 
49152 
40475 
55296 
99967 
110592 
9791 



20480 
5533 
61440 



31766119 

9437184 

43837325 

42467328 

2911646591 

530841600 

820810921 

98304000 

1657839733 

294912000 

21731177 

1 1468800 

693485 



2752512 



296 



It allows us to rewrite eq (4) as 

l)u = u+ — (u - y)((u + vf - 4) - (co - l)Di/ (5) 

8 

From this equation the following parity rules can be read off immediately: 

a) For n even u n has purely real coefficients 

b) For n odd u n has purely imaginary coefficients and <j) n = 
Starting with 

no = 2 £, uj = i(Va f 3 + l /4 £~ 3 ) and o> = 1 

the coefficients co w and the polynomials can be determined recursively by the following method: 
For even n the terms of order e n in eq (5) give rise to the differential equation 

f)u„ = u„ + F tt (£) - 2o> n | 

F n (f ) is a known polynomial since it depends only on //*.(£) and (i) k with A' < n. 

It will contain a linear term with a real coefficient a so that F n ({j) — a + F *(£). Since £ is a solution to the 
homogeneous equation this is a critical term and we have to set co w = a/2. The general solution to the 
remaining equation 

ih,„ = ,/„ + /•'*(£) 

is then 

"« = «£ + "*(f) 

where a*(£) is a polynomial in £. The parameter a will be determined at the next order. The differential 
equation for terms of odd order/? 4- 1 reads 

/;«„h = «■+! + F n+i (i) + «*U + \t 3 - 1' - \t*f- 

F n+ i(£;) denotes again a known polynomial which will contain a linear factor of the form ibij. To avoid secular 
terms we require a = ~b. The remaining equation can then be integrated to give the polynomial U n+1 . The 
solution to the homogeneous solution is set to zero this time out of convenience and to obtain the phase 
shifting solutions of [1]. 

The second author used his algebraic processor POLYPAK to implement the above algorithm by computer 
and to check the results in [1]. POLYPAK is a package of PL/I computer programs for the manipulation of 
real or complex powerseries in several variables. It was derived from MAO. All computations were performed 
on the Amdahl 470/V6 of the University of Cincinnati which runs under the IBM operating system VS2 
release 1.7. 

For the above problem our series had two variables e and £ and the coefficients were stored as complex 
double precision floating point numbers. We determined the solution a and the series a) through terms of 
order 6 4 ° and confirmed the results given in [ 1]. 

Our computations took 57 seconds. This compares with 130 seconds for a FORTRAN program which 
Melvin used to implement his formulas and which he ran on a CYBER 175. A comparison of the two times 
may not mean very much as we compare different machines and different programs. Nevertheless, we believe 
that it indicates that the use of POLYPAK did not generate any unreasonable overhead. The use of an 
algebraic processor is already justified by the ease with which the above algorithm can be coded. 
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